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ABSTRACT 

■ We present the first three-dimensional circulation models for extrasolar gas giant 

■ atmospheres with geometrically and energetically consistent treatments of magnetic 
drag and ohmic dissipation. Atmospheric resistivities are continuously updated and 
calculated directly from the flow structure, strongly coupling the magnetic effects with 
the circulation pattern. We model the hot Jupiters HD 189733b (T cq « 1200 K) and 
HD 209458b (T oq ~ 1500 K) and test planetary magnetic field strengths from to 
30 G. We find that even at B = 3 G the atmospheric structure and circulation of 
HD 209458b are strongly influenced by magnetic effects, while the cooler HD 189733b 
remains largely unaffected, even in the case of B = 30 G and super-solar metallicities. 
Our models of HD 209458b indicate that magnetic effects can substantially slow down 
atmospheric winds, change circulation and temperature patterns, and alter observable 
properties. These models establish that longitudinal and latitudinal hot spot offsets, 
day-night flux contrasts, and planetary radius inflation are interrelated diagnostics of 
the magnetic induction process occurring in the atmospheres of hot Jupiters and other 
similarly forced exoplanets. Most of the ohmic heating occurs high in the atmosphere 
and on the day side of the planet, while the heating at depth is strongly dependent 
on the internal heat flux assumed for the planet, with more heating when the deep 
atmosphere is hot. We compare the ohmic power at depth in our models, and estimates 
of the ohmic dissipation in the bulk interior (from general scaling laws), to evolutionary 
models that constrain the amount of heating necessary to explain the inflated radius of 
HD 209458b. Our results suggest that deep ohmic heating can successfully inflate the 
radius of HD 209458b for planetary magnetic field strengths of B > 3 — 10 G. 
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Introduction 



Hot Jupiters are unlike any planet in our solar sys tem and their atmosph eric circulation exists 
in an entirely new regime (for an extensive review, see lShowman et al.ll2010l ); it has recently been 
recognized that their circulation can also be influenced by the planet's own magnetic field. The 
atmospheres of many of these planets are hot enough that they should be weakly thermally ionized, 
for pressures less than the transition to metallic hydrogen at ~ 1 Mbar and greater than a transition 
to photoionization at ~mbar pressures, with alkali metals providing the primary source of ions. As 
the charged particles in the mostly neutral flow are advected around the planet by very fast winds, 
their circulation through the planetary magnetic field should result in the generation of a secondary 
component to the field an d associated current s. The effects on the atmosphere will be a bulk Lorentz 
force dra g on the winds (Perna et al] 2010a) and localiz ed heating from ohmic dissipation of the 
currents (|Batvgin &: Stevensonll201Cll ; iPerna et al.ll2010bl ). 



There are several ways that magnetic effects may influence observable properties of hot Jupiters. 
One of the most widely recognized is that the ohmic heating may provide the extra source of 



heating required to explain the unexpected la r ge radii of some hot Jupiters (jBatygin et al.l 12011 



Laughlin et al 



2011 



Huang Cumming 



2012 



Wu &: Lithwickll2012l ). For a given planetary mag- 



netic field strength , there should be more heating on planets subject to higher levels of irradiation 
(j Perna et al.1 12012| ) , but these planets will also experience stronger drag on their winds and we 
may expect an anti-correlatio n between the amount of radius inflatio n and the offset of the hot 
spot from the substellar point (|Menoull2012al ; lRauscher &: Menoull2012r). There is evidence for tem- 
perature inversions in the atmospheres of many hot Jupiters (e.g., iMadhusudhan k. Seagerll2010l . 
and references therein) and, while the presence of clouds or s tratospheric absorbers will influence 
the pressure levels at which ohmic heating occurs (lHendl2012l). it is possible that the temperature 
inversions are in fact produced by ohmic heating (|Menou 2012b ). Finally, with future observing 
facilities we may be able to directly measure the upp er atmosphere wind speeds on these pla nets 
and thereby constrain the strength of magnetic drag (jMiller-Ricci Kempton Rauscherll2012l ). 



Due to the complexity of this topic, most of the work so far has estimated ohmic heating rates 
from three-dimensional circulation models or analytic assumptions, without consistently including 
the feedback of magnetic drag and heating on the circulation pattern. We have previously published 
models that attempted to estimate the effect of magnetic drag on hot Jupiter circulation b y in- 



cluding a simple, approximate form for the drag (IPerna et al 



2010a 



Rauscher fc Menou 2012. also 



used for Miller- Ricci Kempton & Rauscher 2012), but here we improve upon that work in several 
important ways. First, instead of estimating the drag strength, we now calculate it directly from 
local conditions (temperature, density) and update it with each timestep, which makes it strongly 
coupled to the dynamics. This also means that the magnetic drag is no longer uniform at each 
pressure level, but instead varies by many orders of magnitude around the planet. In addition, we 
include geometric effects due to an assumed aligned dipole field: only the zonal (east-west) compo- 
nent of the flow experiences drag and the strength of drag is also dependent on latitude (with zero 
drag at the equator). Finally, we convert the kinetic energy lost through drag into heating from 
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ohmic dissipation. By including all of these effects, we present the first atmospheric circulation 
model with geometrically and energetically consistent magnetic drag and ohmic dissipation. 

The obvious unknown in studying this topic is the strength of hot Jupiter magnetic fields. 
We can use scaling laws to predict field strengths, but our knowledge is limited by the unknown 
complexity of planet interiors and an incomplete understanding of dynamo theory (see reviews by 
Christensenll2010l ; IStevensonll2010l ). Nevertheless, such scaling la ws estimate hot Jupiter mag netic 
field strengths to be anywhere from a few to tens of Gauss (see lReiners &; Christensenll2010l . and 
references therein) . Of course it would be preferable to actually measure the magnetic field strength 
for any particular planet, but unfortunately this will most likely require the use of indirect methods 
in which the signature of the planet is buried in the stellar signal. These include the measurement 
of changes in the stellar chromospheric or X-ray emission due to interaction between its and an 
orbiting planet's magnetic field (e.g., Shkolnik et al. 2005, 2008, Kashyap et al. 2008; although 
see Poppenhaeger et al. 2010, 2011, Miller et al. 2012), dire ct detection of radio em ission from the 
planet (which is likely too weak for current capabilities, see iGrieBmeier et al.ll201ll . and references 



therein), and the signature of a planet's magnetic field in its influence on the rate and geometry of 



Yelk 


2004; 


Adams 


2011; 



Trammell et al. 



2011 



We describe our numerical model in Section [2j with particular attention to the new imple- 
mentation of magnetic effects. In Section [2.21 we catalog the set of models presented in this paper, 
explaining our choices for the range of physical parameters used. We begin the analysis of our 
results with a detailed look at the model of HD 209458b with B = 3 G (|3.ip . followed by an 
examination of models with increasing magnetic field strengths (j3.2|) . and a comparison between 
our models of HD 209458b and the cooler planet HD 189733b (|3.3j) . We then discuss the influence 
of numerical resolution on our results (|3.4p . In Section 13.51 we examine the observable properties of 
our models and in Section [3.61 we discuss global ohmic heating rates, commenting upon how they 
could influence a planet's thermal evolution and its radius. We conclude with a summary of our 
main findings in Section [H 



2. Our Numerical Model 



We have updated our General Circulation Model (GCM) to include the effects of geometrically 
and energetically consistent magnetic drag and heating. Our G CM is adapted from t h e Inte rmediate 
General Circulation Model (IGCM) originally developed by iHoskins &; Simmons! (|1975l ). with a 
pseudo-spectral dynamical core that solves the primitive equations of meteorology. We disc u ss ou r 
ada ptation of the code t o stud y synchronously rotating gas giants in iMenou k, Rauscherl ([2009) 



and iRauscher &; Menoul (120101). and our implementa tion of a standard two-stream, double-gray 
radiative transfer scheme in IRauscher &; Menoul (|2012l ). hereafter RM12. For our radiative transfer 
scheme, we separate the optical and infrared wavelengths, with the attenuation of the incident 
stellar (optical) flux controlled by a constant optical absorption coefficient, while the absorption 
and emission of thermal radiation at each level is set by a constant infrared absorption coefficient. 
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2.1. Modeling the magnetic effects 

We employ several simplifying assumptions in order to mode l complex magnetic p rocesses with- 



out using a full magnetohydro dynamics (MHD) calculation. In iPerna et al.l (|2010al ) we compared 
the relative importance of the non-ideal MHD terms in the induction equation to determine that 
most of the atmosphere is within the purely resistive MHD regime and that magnetic Reynolds num- 
bers are generally much less than 1, although there are localized regions where these assumptions 
can break down. We additionally assume that the planet's magnetic field (presumably generated 
deep within the interior) has a dipole geometry and is aligned with the planet's rotation axis; it 
also remains unaltered by the weakly ionized flow in the atmosphere. 

We only apply magnetic drag to the zonal (east-west) wind (u), leaving the meridional (north- 
south) wind unaltered, by adding a term du/dt = —u/T ma _ g to the momentum equation. In the case 
of an aligned dipole field, any drag on the meridional flow is not significant until near the poles, 
where the field becomes more radial. As a simplification we do not vary the magnetic field strength 
along the planet surface — as it should for a true dipole — effectively ignoring the magnetic geometry 
at the poles. We choose to use a constant field strength instead of more complex geometry in part 
because the detailed geometry of the actual planetary field is unknown. However, the meridional 
flow can be fairly significant over the poles, especially high in the atmosphere. Future work will 
be required to expand the magnetic formalism we use and to determine how large of an impact 
meridional drag could have on the circulation. 

The magnetic timescale (T mag ) is calculated from the chosen magnetic field strength (B) and 
local conditions (density, p; temperature, T; and latitude (ft): 

r mag (B,p,T^) = W j^ ] ■ (1) 

(For a derivation of r mag , see Perna et al. 2010a.) Note the geometric dependence of r mag ; due to 
the assumptions of latitudinal currents and an aligned dipole field, jx B ' = at the equ ator (eft = 0) 



and there is no drag on the flow. The resistivity (77) is calculated as in iMenoul (|2012al ): 

7] = 230\/T/x e cm 2 s" 1 (2) 

with the ionization fraction {x e <C 1) calculated from a form of the Saha equation that takes into 
account the first ionization potential (e^) of all elements from hydrogen to nickeljl assuming solar 
abundance for each element (a^): 

ni = ( — ) n (3) 



X 2 - 1 /"W, \ 3 / 2 



x 2 ■ riikT \ h 2 



( ^ (W 2 



e.i 



(kTf^expi-eJkT) (4) 



2 By using the first 28 elements, rather than just potassium (as in Perna et alj|2010a ) , the resistivities and timescales 
are decreased by a factor of ~2, for the conditions of interest here. 
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i=l,28 



(5) 



where n is the number density (= p/p, with p being the mean molecular weight). The local 
resistivities and magnetic timescales throughout the entire atmosphere are updated at each timestep 
in the simulation, making the magnetic effects strongly coupled with the atmospheric circulation. 

We assume that all of the kinetic energy lost via magnetic drag is returned to the atmosphere 



(Liu et al. 


2008; 


Perna et al. 


2010a: 


Menou 


2012a) 



1 4vrr/ , 2 
.2 3 



,5 2 



sin i 



u 



dt J ohm P 61 " A7T VP 'mag 

which is energetically consistent with our form for the drag (du/dt = —u/r n 



(6) 



In Figure Q] we plot the electrical resistivity (77) as a function of temperature and density, 
over the range of values expected for hot Jupiter atmospheres. Although only weakly dependent 
on density, the resistivity is a strong function of temperature. A second plot in the same figure 
shows an example of how the magnetic timescale (T mag ) varies throughout a planet's atmosphere. 
Timescales are generally shorter deeper in the atmosphere, but the strongest variation — by many 
orders of magnitude — is in the upper levels, between the hot dayside and the cold nightside. FigureQ] 
also shows the geometric dependence of r mag , resulting in differing timescales at the same values of 
the local density and temperature. 

It is impor tant to note that we are applying a formalism derived assuming axisymmetry 
(|Liu et al.ll2008l ) to a non-axisymmetric problem. Figure [1] shows the huge difference in resistivities 
between the day and night sides of the planet. In addition, there are regions of the atmosphere 
(usually at low pressure) where the circulation is not dominated by a zonal jet (or jets), meaning 
that there is a significant non-axisymmetric component. As we will see, this can become even 
more of an issue once the magnetic effects are included in the models (for example, see the flow 
patterns at low pressure in Figures [2] and H]) . While the use of this formalism may be justified in 
first approximation, we would (as always) benefit from a fuller, non-axisymmetric theory. 

In order to maintain numerical stability we prevent extremely strong drag or heating by setting 
a minimum magnetic timescale, T magjm i n , which is used instead of the timescale calculated from 
Equation [1] if it is greater than r mag . In the models presented here we set the minimum magnetic 
timescale equal to our hyperdissipatiorjf] timescale (= 0.005 in units of the planet's rotation period, 
P ro t). For comparison, the timestep used in our models is ~ 0.0002 P ro t- Among the range of 
parameters we tested, it was only the model with strongest magnetic field strength (B = 30 G) 
and the hottest atmosphere (HD 209458b) where our choice of T magjimn = 0.005 limited the full 
strength of the magnetic effects (and that model failed to run to completion, see Section [3.21) . 



3 Hyperdissipation is a common technique used to prevent the build up of numerical noise on the smallest resolved 
scales. 
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Fig. 1. — Left: the electrical resistivity of the atmosphere, 77, as a function of the local temperature 
and density, over the range of values expected in hot Jupiter atmospheres. The contour levels give 
log 10 (??) in cm 2 s _1 . Right: the timescale for magnetic drag and heating, r mag , for our model of 
HD 209458b with B = 3 G. The color of each point gives log 10 (r mag ) in seconds. Each group 
of points (roughly vertically aligned) corresponds to one of the discrete pressure levels in the 
simulation. The variations in r mag for a given density and temperature are due to the geometric 
dependence of r mag (see Equation [T]) . 



We initialized each model with the atmosphere at rest (no winds). The initial temperature 
structu r e was set to be uniform at each pressure level, using a temperature-pressure profile from 



Guilloti (120101 ) with parameters matching those used for our radiative transfer, and with the av- 



eraging parameter set to / = 0.375, between a global average (/ = 0.25) and a dayside average 
(/ = 0.5), in order to minimize the time needed for the atmosphere to equilibrate. (A global 
average was not chosen because it would have led to very strong heating in the substellar regions.) 
The atmosphere was then allowed to heat and cool according to our radiative transfer scheme, 
accompanied by the acceleration of winds. After 10 orbital periods we introduced the magnetic 
effects of drag and heating, linearly increasing from zero to reach their full values at 20 orbital 
periods. We tested using 2/4 P rot or 50/100 P rot instead of 10/20 P rot for the implementation of 
magnetic effects. By waiting to apply the magnetic effects until later in the run, the winds were 
able to accelerate to faster values early in the simulation, but by 1000 P ro t we found no significant 
difference between the models. 



2.2. Our set of models 



We have tested the effects of magnetic drag and heating over a range of physical and numerical 
parameters. The default resolution we used was T31L30, the horizontal spectral resolution corre- 
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sponding to ~4° and the 30 vertical levels logarithmically spaced in pressure, from 1 mbar to 100 
bar. These are the same parameters we used in the models without magnetic effects from RM12 
and we have also chosen the same hyperdissipation strength (V 8 , r = 0.005 P TO t) and the same 
length for the runs (out to 2000 P ro t= -forb)- In Section [3~41 we will discuss our tests of numerical 
resolution. 

We ran models for planetary magnetic field strengths of B = 0, 3, 10, and 30 G. As discussed 
in Section [H it is not clear what field strength to expect for Jupiter-mass planets on orbital periods 
of a few days, so our values span a range from weak to strong fields. We used planetary param- 
eters meant to represent two well known hot Jupiters: HD 189733b and HD 209458b, as shown 
in Table [TJ In addition to being the two hot Jupiters with the most measurements, the difference 
in equilibrium temperature between these two planets (AT cq ~ 300 K) means that we can test 
two levels of thermal ionization, with a difference of ~2 orders of magnitude in electrical resistivi- 
ties. Although the physics of radiative transfer is highly detailed and differences between planets, 
particularly in composition, can lead to variation in their radiative properties, here we choose a 
more straightforward approach and use the same optical and infrared absorption coefficients in for 
both planets. The hotter atmosphere of our HD 209458b model is solely the result of a higher 
incident stellar flux than for the HD 189733b model. Note that we use different planetary radii, 
surface gravities, and rotation rates for these two planets, as appropriate, but these differences are 
secondary compared to the disparate levels of stellar heating. 

The level of thermal ionization in a planet's atmosphere is dependent both on temperature 
and on the abundance of elements with low ionization potentials (particularly alkali metals). Since 
we may expect planets to be enriched in metals (e.g. Jupiter has abundances ~ 3x solar, Wong 
et al. 2004), this is an additional physical parameter that we must consider in our study of 
magnetic effects. The uncertainty in atmospheric composition factors into the resistivity and may 
be somewhat balanced by our uncertainty in the magnetic field strength, since these parameters 
together control the strength of the magnetic effects (Equation [I]) . Since we expect magnetic effects 
to matter less in the atmosphere of HD 189733b than the hotter HD 209458b, we tested whether we 
could bring them to the same level by running a model of HD 189733b with B = 30 G and metal 
abundances increased 3x above solar, a value chosen somewhat arbitrarily, in lieu of a consensus 
on the measured metallicity for this planet. 

The physical effects that result from increased atmospheric metallicity are: 1) an increased 
abundance of ions from thermal ionization, 2) an increase in the mean molecular weight (MMW), 
and 3) a change in opacities throughout the atmosphere. While the slight increase in MMW could 
have a small effect on the atmospheric dynamics, changing the opacities is known to have a large 



effect (|Dobbs-Dixon &: Linl 120081 : IShowman et al.l 120091 : iLewis et al.l l2010l ) , since this controls the 



structure of radiative heating throughout the atmosphere. However, here we are only going to adjust 
the elemental abundances in the ionization calculation, without changing the other parameters of 
the model, in order to isolate the effect of increased ionization. Since we are using two constant 
absorption coefficients for our radiative transfer, and these were chosen in order to roughly match 
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Table 1. Model parameters used 



Parameter HD 189733b HD 209458b Units 



Radius of the planet, R p 


8 x 10 7 


1 x 10 s 


m 


Gravitational acceleration, g 


22 


8 


m s 


Rotation rate, Q, 


3.3 x 10" 5 


2.1 x 10~ 5 


s- 1 


Corresponding period, P rot = 2tt/Q 


2.2 


3.3 


day© 


Equilibrium temperature 21 , T eq 


1200 


1500 


K 


Incident flux at substellar point, i^vis^rr 


4.74 x 10 5 


1.06 x 10 6 


W in" 2 


Corresponding temperature, Tl rr 


1700 


2078 


K 


Internal heat flux, -FfiR^nt 


3500 




W m- 2 


Corresponding temperature, T; nt 


500 




K 


Optical absorption coefficient, k v i s 


4 x 10" 


-3 


2 —1 

cm g 


Infrared absorption coefficient, KrR 5 o 


1 x 10" 


-2 


2 —1 

cm g 


Infrared absorption powerlaw index, a 









Specific gas constant, TZ 


3523 




J kg" 1 K- 1 


Ratio of gas constant to heat capacity, TZ/ cp 


0.286 







a assuming an albedo of zero 
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expected temperature-pressure profiles from ID models, there is no clear prescription for how we 
would change our opacities to match the change in metallicity. An interesting avenue for future 
work will be atmospheric circulation models that investigate the dual influence of metallicity on 
radiative heating and on magnetic effects. 

Finally, we test the importance of the internal heat flux (characterized by Tj nt ) on the amount 
of ohmic heating in the atmosphere. Our choice for Tj n t does not affect the atmospheric circulation 
high in the atmosphere, near the photosphere (as demonstrated for non-magnetic models in RM12). 
However, whether ohmic heating could lead to radius inflation depends not just on the strength of 
heating, but also on the depth at which it occurs, with d eeper heating more likely to influence the 



global structure of the planet (jGuillot Showmanll2002l ). The deep atmosphere will be hotter for 



higher values of T; nt and we expect that this should result in more ohmic heating. Although a full 
study could be made examining a range of Tint, we will leave that for future work and here simply 
compare B = 3 G models of HD 209458b with T int = 500 K (our default) and T int = 100 K. 



3. Results 

We begin the presentation of our results with a detailed look at the model of HD 209458b 
with B = 3 G in order to demonstrate how magnetic effects function in hot Jupiter atmospheric 
circulation. Aside from the inclusion of magnetic drag and heating, this model is identical in all 
parameters to one presented in RM12, allowing for a clean comparison against a non-magnetic 
model. 



3.1. HD 209458b with B = 3 G 

The first difference between this model and the non-magnetic version in RM12 is in the develop- 
ment of the flow, starting from an initial condition with no winds. In the absence of magnetic drag 
the atmosphere accelerates, first in the directly forced region of the atmosphere (above ~1 bar) and 
then, through the vertical transfer of momentum, in the levels below the photosphere. Eventually 
the winds reach their peak speeds and the atmosphere is in a quasi-steady state. We find that the 
presence of magnetic drag alters this pattern. While there is still an initial acceleration of winds 
in the upper and lower levels, this is quickly followed by the response of the magnetic drag, decel- 
erating the winds. This initial ramp-up/ramp-down period lasts ~1200 P ro t for the HD 209458b 
B = 3 G model, after which follows repeated acceleration and deceleration of the winds, with no 
strict periodicity, but a rough timescale of ~100 P ro t (Figure [8] shows this behavior, in the zonally 
averaged flow at the equator). We find that the entire atmosphere gains and loses kinetic energy 
(and total enthalpy) with these variations, although by the end of the 2000 P ro t simulation the 
variations in total energies are less than ~1%. 

In addition to slowing the average wind speeds by ~1 km s -1 (see Figure [TT|) . the presence 
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of magnetic drag also disrupts the development of the strong super-rotating equatorial jet seen 
throughout the atmospheres of most hot Jupiter models (se e Figure [3j), by inhibiting the jet- 



pumping mechanism identified by IShowman Polvanil (|201ll ) . This has important consequences 



for the temperature structure of the atmosphere, as seen in Figure [2j where we show horizontal 
slices at different depths in the atmosphere. (We also plot the specific ohmic heating rates in this 
figure, calculated from the local conditions as per Equations [T] through [6j) 

The temperature structure in this model is more strongly tied to a hot day/cold night pattern 
than it is in the non-magnetic model. Whereas in the non-magnetic model day-night temperature 
differences decrease with increasing pressure, and are negligible below ~1 bar (see Figure 2 of 
RM12), in this magnetic model temperature contrasts of 1000 K are maintained down to at least 
1 bar, only becoming homogenized below ~4 bar. The simplest explanation for this is that the 
slower wind speeds mean that gas heated on the day side has more time to cool before reaching 
the night side. This is valid in the highest regions of the atmosphere, where we see the same type 
of substellar-to-antistellar (SSAS) flow as in the non-magnetic model; however, slightly different 
explanations are required at deeper levels. 

If we compare the flow at the infrared photospheres of this and the non-magnetic model (see 
Figure 5 of RM12), we can see that in this magnetic model the hottest region of the atmosphere 
remains at the substellar point, instead of being advected eastward by tens of degrees in longitude. 
At these pressures the difference in temperature structures is not just the result of slower winds, but 
also the lack of a strong eastward equatorial jet. Even though there is no drag along the equator 
(see Equation [TJ, the super-rotating jet is inhibited because the magnetic drag acts on the zonal 
(but not meridional) flow, and p revents the slanted, up-grad ient transport of angular momentum 



toward the equator described by IShowman &: Polvanil (|201ll ). Instead, we find a flow that is still 



mostly SSAS with a convergent feature at ~135° E, as at higher levels. 

This convergent featurejll seen across multiple levels, is a feature in common with the non- 
magnetic model. As in that model, here we find that downward motion due to the convergent flow 
results in adiabatic heating at that location in deeper pressure levels. In the non-magnetic model 
this local hot spot is superimposed on the hot regions advected eastward from the day side. In this 
model, however, the hot spot remains an independent feature and is in fact enhanced by ohmic 
heating, especially at layers below the optical photosphere (~1 bar), where stellar heating is no 
longer important. Although this leads to an effective eastward shift in the hottest regions of the 
atmosphere, note that this is a distinctly different phenomenon from eastward advection due to an 
equatorial jet. 

Finally, although the 135° hot spot is an exception, the temperature structure in this model is 
influenced more strongly by the presence of magnetic drag, than by the (coupled) effect of ohmic 



4 See Section 4.3 of Rauscher fc Menoul for a discussion of this feature, found in many circulation models, 



including those from other groups (e.g.. IShowman et al.ll2009l : iDobbs-Dixon et al.ll2012l ) 
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Fig. 2. — Ohmic heating (shown as blue contours at 20, 40, 60, and 80% of the peak value, and 
coincident with the magnetic drag) for the model of HD 209458b with B = 3G, plotted with winds 
(arrows) and temperatures (color scale, in K). These horizontal slices through the atmosphere (in 
cylindrical projection, centered on the substellar point) are at pressure levels of 8 mbar (top left), 56 
mbar (near the infrared photosphere, top right), 260 mbar (middle left), 800 mbar (middle right), 2 
bar (bottom left), and 3 bar (bottom right), with peak wind speeds of 8.1, 5.9, 4.6, 2.5, 1.2, and 1.3 
km s , and peak ohmic heating rates of 22, 3.4, 1.6, 0.10, 0.014, and 0.0039 W kg -1 , respectively. 
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heating. When we compare rates of magnetic and radiative heating throughout the atmosphere, we 
find that almost everywhere the radiative heating dominates. Deep in the atmosphere the ohmic 
heating becomes a non-negligible fraction of the radiative heating, but we will delay discussion of 
this point until our section about global heating rates (|3,6p . 



3.2. Changes in the circulation due to increasing magnetic field strength 

We do not know the strength of hot Jupiter magnetic fields. In our previous work on this 
topic we considered models of HD 209458b with magnetic field strengths of B = 3, 10, and 30 G 
rtPerna et al.ll2010al lbl. :IM12) and we have used those same values here in order to facilitate com- 



parison. We have improved upon that previous work by including ohmic heating in addition to 
drag, by applying drag only to the zonal component of the wind, and by continuously calculating 
the strength of the magnetic effects based on changing local conditions, rather than just using a 
single value of r mag for each pressure level. Using this updated code, we find that we are unable 
to successfully run models of HD 209458b with B > 20 G, for reasons we will discuss below. To 
begin, we will compare our models of HD 209458b with B = 0, 3, and 10 G. 

I n our previous mod els we saw a gradual change in circulation as B was increased from to 



30 G (jPerna et alj|2010al . RM12), but now with our more realistic scheme we see a sharp change in 
the circulation from the no-drag to the B = 3G model and then a more subtle change as the field 
strength is increased from B = 3 to 10 G, as shown in the zonal wind profiles plotted in Figure [3] 
The B = 10 G model has similar circulation patterns as the 6 = 3G model, although with slightly 
slower wind speeds and with more of a departure from hemispheric symmetry (this can also be 
seen in Figure ITT]) . primarily due to an increased perturbation of the equatorial jet, causing it to 
meander farther to the north and south. 

These same trends continue to higher magnetic field strengths, as seen in the developing flow 
of the B = 30 G model, before it becomes numerically unstable and crashes. Even with 1.5 x more 
timesteps per P ro t> compared to the other models, this model would only run to 271 P rot before 
crashing. (We also tested using B = 20 G instead, but that model was likewise unable to run for 
the full 2000 Prot-) Although this uncompleted run cannot be used for a reliable picture of the 
circulation at B = 30 G, it is nevertheless informative to examine the atmosphere before the crash 
in order to understand the limitations of the model. 

In Figure H we show layers of the HD 209458b B = 30 G model, at 271 P rot . At the 3 mbar 
level we see the effect of artificially limiting the minimum value for r mag in our models. Here the 
ohmic heating should be a combination of that from the zonal flow across the terminator (where 
winds are fast and r mag > r magim i n ) and discrete regions on the day side (where zonal winds are 
almost nonexistent, but r mag is very short, < T magimrn ). However, the lack of this heating is only a 
small loss; integrating the heating rates over this pressure level (see Equation IA2j) . we find that it 
would have added only 2 x 10 18 W to the total ohmic heating in this pressure level, = 6 x 10 19 W. 
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Fig. 3. — Plots showing the zonal average of the zonal (east-west) wind (in m s _1 ), as a function of 
latitude and pressure, for our models of HD 209458b with B = 0G (left), B = 3 G (middle), and 
B = 10 G (right). The yellow line separates eastward (positive) from westward (negative) flow. 
These are from snapshots of each model at 2000 P ro t- 
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Fig. 4.— Snapshots of the HD 209458b B = 30 G model at 271 P rot , immediately before the 
simulation became numerically unstable and crashed. The dark blue contours show the ohmic 
heating (coincident with the magnetic drag), the light blue contours show the ohmic heating that 
would have occurred if we had not set a limit on its strength (r magim i n = 0.005), the color scale gives 
the temperature (in K) and winds are represented as arrows. Contour levels are plotted at 20, 40, 
60, and 80% of the peak value. Left: 3 mbar pressure level, with a peak wind speed of 8.1 km s _1 , 
a peak ohmic heating rate of 1000 W kg -1 , and a peak in "missing" heating of 52 W kg -1 . Right: 
380 mbar level, with a peak wind speed of 3.1 km s _1 and a peak heating rate of 43 W kg -1 . 

The plot on the right side of Figure H] shows a slice of the atmosphere at 380 mbar. Here we 
see a significant break in hemispheric symmetry, with a perturbed equatorial jet that only exists on 
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the night side and is strongly dragged as it tries to flow around to the day side. In fact, the ohmic 
heating that results from this drag has a peak value of 43 W kg -1 , which is greater than 10% of 
the local rate of heating from the stellar insolation (as calculated by our radiative transfer scheme). 
Below the optical photosphere, at 2 bar, there is also a region of strong ohmic heating associated 
with the hot feature at 135° and there the local ohmic heating is 2.8 W kg -1 , or ~2% of the net 
radiative cooling. Since our calculation of resistivities and magnetic timescales are updated with 
each timestep, the ohmic heating is strongly coupled to the dynamics and these localized regions 
of intense heating will fluctuate with changes in the flow, in contrast to the more steady stellar 
heating. These intense local heating rates may be the cause of our numerical instability and point 
to the drastic influence of a B = 30 G field on the atmospheric circulation of a planet as hot as 
HD 209458b. 



3.3. Comparison between HD 189733b and HD 209458b 

The main distinguishing feature between HD 189733b and HD 209458b, at least as applies 
to our models here, is the ~300 K difference between their equilibrium temperatures and the 
resulting 2 orders of magnitude difference in their atmospheric resistivities. This alone indicates 
that magnetic effects should be far less important for HD 189733b than we found them to be for 
HD 209458b. 

Since we have not previously published any non-magnetic models of HD 189733b, we first briefly 
present one here. Figure [5] shows several horizontal slices through the atmosphere, demonstrating 
that the circulation on HD 189733b is qualitatively similar to HD 209458b: there is a strong 
eastward equatorial jet that extends throughout most of the atmosphere (the zonal wind profiles 
are very similar, except with slower peak values for HD 189733b) and the jet becomes increasingly 
efficient at advecting the hottest regions away from the substellar point at deeper pressure levels. 
At the photosphere this results in an eastward shift of the hot spot by ~20°, and by pressures of a 
few bars the temperatures are well homogenized around the planet. Although the winds are slower 
on HD 189733b than on HD 209458b, the cooler temperatures mean that the radiative timescales 
are also longer, and HD 189733b has lower temperature contrasts than HD 209458b as a result. 

A comparison between models of HD 189733b with B = 0, 3, 10, and 30 G shows very minimal 
differences. Wind speeds are slightly slower for higher B models, but the strong equatorial jet 
remains and the general structure of the atmosphere is largely unchanged. In Figure [6] we plot 
the photospheres of models with B = 30 G, one where we have assumed solar metallicity and one 
where we have used 3x solar. Even in the case of 3x solar metallicity the atmosphere looks almost 
identical to the B = model. The increase in metallicity by a factor of 3 leads to a decrease 
in atmospheric resistivities by ~3, which is almost equivalent to an increase in the magnetic field 
strength by a factor of y/3 (see Equation [[]) . Our results show that magnetic effects seem unable to 
significantly alter the atmosphere of HD 189733b, even at super-solar metallicities and for strong 
magnetic field strengths. 
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Fig. 5. — Temperature and wind maps for our HD 189733b model with B = G. The pressure levels 
shown are 8 mbar (top left), 180 mbar (the infrared photosphere, top right), 1 bar (bottom left), 
and 6 bar (bottom right), with maximum wind speeds of 7.0, 6.1, 5.2, and 4.1 km s -1 , respectively. 
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Fig. 6. — Temperature, wind, and ohmic heating maps at 180 mbar for our HD 189733b model 
with B = 30 G and solar metallicity (left) or 3x solar (right). The maximum wind speeds are 5.6 
and 5.3 km s _1 and the peaks in ohmic heating are 5.9 and 8.4 W kg -1 , respectively. 
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3.4. Tests of numerical resolution 

We might expect that the magnetic effects we have included are especially sensitive to res- 
olution; their strength (via r mag , Equation [T|) is dependent both on the latitude (sampled more 
frequently at higher resolution) and on the detailed temperature structure, with a small contrast in 
temperature resulting in an exponential difference in r mag . The ohmic heating and drag patterns 
in Figure [2] show both of these factors at work. Often the magnetic effects are strongest near the 
equator, due to the higher temperatures and stronger winds commonly found there, but careful 
examination shows that the effects decrease sharply at the exact equator. It is therefore beneficial 
to test whether the horizontal resolution of the simulation has a strong influence on the results. 

However, a comparison of models at different resolutions becomes quickly complicated by 
the effect of hyperdissipation. This common numerical scheme is used to reduce the build-up of 
energy at the smallest resolved scales by applying a high-order operator to the divergence, relative 
vorticity, and temperature fieldsll Hyperdissipation is meant to represent subgrid processes that 
will continue the cascade of energy down to smaller scales, where kinetic energy is eventually 
dissipated. Unfortunately it is not possible to calculate the appropriate strength and order for 
hyperdissipation from physica l principles alone; ideally it should be carefull y tested for application 



to each particular model (e.g. iThrastarson &: Choll201ll ). iHeng et al.l (|201ll ) demonstrated that for 



hot Jupiter models at a given resolution, the maximum wind speeds were dependent on the choice 
of hyperdissipation strength, and, similarly, that maximum wind speeds vary when an identical 
hyperdissipation strength is applied to models at a range of resolutions. Since magnetic effects 
are sensitive both to the speed of atmospheric winds and to the spatial structure of temperature 
and wind patterns, this makes the relationship between resolution, hyperdissipation, and magnetic 
effects a difficult problem to untangle. 

We ran tests at lower horizontal resolution, T21 (~5.6°), and a couple of limited cases at 
T42 (~2.8°, and computationally very expensive). Although the best choice for hyperdissipation 
strength should generally change with resolution, it could also depend on the magnetic effects and 
how they influenced the circulation pattern. In lieu of performing a comprehensive resolution study 
we kept the hyperdissipation identical to that used for our T31 models. From our resolution tests 
in RM12, we know that this did not result in too much over- or under-damping in our drag-free 
models at T21 and T42. 

In general we found no significant differences between equivalent models at different resolutions; 
the resulting circulation patterns, temperature structures, and observable properties were fairly 
similar. In Figure [7] we compare results from a sample of our tests, both for models in which 
magnetic effects were nonexistent or had a weak effect on the circulation (HD 209458b with B = G, 
and HD 189733b with B = 30 G and 3x solar metallicity) , and for models where the magnetic 



5 See RM12 for a description of how we chose the hyperdissipation to use (raiss = 0.005 P rot and V 8 ) for our 
drag-free model of HD 209458b. 
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effects strongly influenced the flow (HD 209458b with B = 3G, and with B = 10 G). As expected, 
the models with no/weak magnetic effects have slower wind speeds in the T21 versions than at T31, 
due to the same hyperdissipation strength being used for both. For the HD 189733b model this 
results in decreased ohmic heating near the equator, since the jet has less kinetic energy available 
to be dissipated through drag and converted to heat. The HD 209458b B = 3 G T21 model 
also has less near-equatorial heating than its T31 counterpart, although the difference is not as 
much as between the HD 189733b models. At B = 10 G the disparity between the T21 and T31 
models of HD 209458b is comparable to the temporal variation within the T31 model. These 
results indicate that as the influence of magnetic drag on the circulation increases, the influence of 
hyperdissipation decreases. It also demonstrates that it is the interdependent relationship between 
resolution, hyperdissipation, and magnetic drag — and not each effect on its own — that influences 
the atmospheric circulation in these models. While this is not an exhaustive test of numerical 
resolution, these results indicate that there is no sudden change in the circulation pattern at higher 
or lower resolution, for a range of physical properties, supporting our use of T31 for the models 
presented in this paper. 

Although the end results were similar, we did notice an interesting peculiarity in the initial 
evolution of the T21 version of our HD 209458b B = 3 G model, compared to the T31 version. 
In Figure [8] we plot the behavior of the zonal mean flow at the equator for models of HD 209458b 
with B = 3 G, over all pressure levels, as function of time in the simulation. The T21 and T31 
versions both initially develop strong eastward equatorial flow throughout the atmosphere, which is 
subsequently decelerated by the magnetic drag. In the T21 model this deceleration happens more 
quickly and results in a flip in the mean equatorial flow, from eastward to westward. (This flip in the 
flow direction also occurs at higher latitudes, so that the dominant global flow becomes westward.) 
The magnetic drag then works to oppose the westward flow and eventually the circulation returns 
to being primarily eastward at the equator, followed by the same irregularly periodic acceleration 
and deceleration of the flow seen in the T31 model. The only other model from our set that exhibits 
this same flip between eastward and westward flow is our T31 model of HD 209458b with B = 3G 
and Tint = 100 K (also shown in Figure [8]). We emphasize that the final circulation states of these 
models are all similar, but we present this peculiarity in order to demonstrate that the presence of 
magnetic effects can influence not just the state of atmospheric circulation, but also its development 
during a simulation. The inclusion of magnetic effects expands the dynamical possibilities for the 
atmospheric flow and introduces greater variability in its development. The potential richness of 
expanded dynamical regimes will be explored in future work. 
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Fig. 7. — Profiles of wind speeds and ohmic heating rates, as a function of latitude, for models 
at T21 and T31 horizontal resolutions. The wind speeds are calculated as the root-mean-square, 
and the heating rates as the sum, over all longitudes and pressure levels for a given latitude. (See 
Appendix A for the calculation of the heating rates.) The profiles are for the models of HD 209458b 
with B = 0G (black), B = 3 G (red), B = 10 G (blue), and of HD 189733b with B = 30 G and 
3x solar metallicity (purple). For each model we show profiles at 2000 P ro t for T21 (dashed) and 
T31 (solid), with an additional profile from one rotation period earlier in the T31 runs (dotted). 
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Fig. 8. — Zonally averaged zonal wind speed (in m s _1 ) at the equator, as a function of depth in 
the atmosphere and time in the simulation (where a "Planet Day" is one rotation period, P TO t)- 
The white line separates eastward (positive) from westward (negative) values. The models shown 
are all for HD 209458b with B = 3 G, with a horizontal resolution of T21 (top), T31 (middle), 
and T31 with T; n t = 100 K (bottom). The values at the equator are calculated as the average of 
±1.86° for T31 resolution and ±2.77° for T21 resolution. 
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3.5. Effects on observable properties 

Most of the observable properties of an exoplanet can contain signatures of its atmospheric cir- 
culation. As the planet transits the star, the detailed shape of the transit curve will depend on the 



geometric shape of the planet, which can be altered by the atmospheric circulation (jBarnes et al 



2009; iDobbs-Dixon et al.ll2012l ). The depth of the transit as a function of wavelength will depend 
on the composition and temperature profiles along the terminator, a region particularly vulnerable 



to the influence of wind s blowing from day to night ([Burrows et al. 



2010 



Fortnev et al.1 [2010b; 



Dobbs-Dixon et al.ll2012l ). With future observational facilities it may be possible to directly mea- 
sure the speed of t h ese winds, as they produce a Doppler shift in spectra taken during transit 
(jSnellen et al.1 hold ; iMiller-Ricci Kempton Rauscherl bold ; Ishowman et al] boij ). The atmo- 
spheric circulation can even have an indirect effect on the size of the planet, via heating by ohmic 
dissipation, as we discuss in the next section. 

In addition to observing the shadow of the planet (and its atmosphere) during transit, we 
can also measure the light emitted by the planet, which for hot J up iters typically peaks in the 
infrared. This is a property of the planet that is self consistently predicted by our numerical 
code. In Figure we plot maps of the infrared flux emitted from the top of the atmospherejf] for 
models of HD 189733b and HD 209458b with B from to 30 G. It is immediately apparent from 
these maps that magnetic effects significantly change the observable properties of HD 209458b, but 
not necessarily those of HD 189733b. We integrate these maps over the planet disk, for viewing 
orientations along the equator, to calculate the change in infrared emission that would be measured 
by a distant observer as the planet orbits its star, assuming that the planet's orbital and rotation 
periods are the same (as we have in our model set-up). These orbital phase curves are plotted in 
Figure [10| where we also plot the phase curve for each magnetic model divided by the curve for 
the 6 = G model of that planet, in order to see small differences between similar curves. 

As expected from the maps in Figure [9j we see a significant change between the B = G and 
B > G models of HD 209458b. The B = model has an effective longitude of peak emission 
that is shifted away from the substellar point (by #Fmax = 12°, as reported in Table [2]), while the 
B > models all have their peak emission well aligned with the substellar point (#Fmax = — 3°). 
There is a lower ratio between the minimum and maximum flux for the B > models than the 
B = one (F m i n / F ma _ x = 12 — 13%, compared to 17%) and the differences between these models 
are greater than the temporal variation in emission within any of the models, which is at a level of 
less than 1%. On the other hand, there is practically no difference between the phase curves for our 
various models of HD 189733b. They all have almost the same effective longitude of peak emission 
(#Fmax = 17°) and flux ratio between minimum and maximum emission {F m \ n / F ma _ x = 21%). The 
relative differences between the phase curves of these models is comparable to the low level of 
temporal variation in emission within each model (<2%). 



6 Note that in our modeling scheme all of the thermal emission of the planet is captured in a single band, so that 
we have no wavelength information. 
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Fig. 9. — Cylindrical maps of the infrared flux (in W m 2 ) emitted from the top boundary of the 
models for HD 189733b (left column) and HD 209458b (right column) with B = G (top row), 
B = 3 G (second row), B = 10 G (third row), and B = 30 G (bottom row). The substellar point 
is at the center of each plot. 
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Fig. 10. — Emitted infrared phase curves for models of HD 189733b (black) and HD 209458b (red) 
with B = G (solid), B = 3 G (dotted), B = 10 G (dashed), and B = 30 G (dash-dotted). 
The bottom panel gives the emitted flux, averaged over the hemisphere facing the observer, while 
the top panel gives the ratio between each magnetic model and the B = G model. The planet 
transits the star at a phase of and passes behind the star at a phase of 0.5; the dip in light during 
secondary eclipse (not shown) provides a measure of the emission from the day side of the planet. 
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Even though HD 209458b is one of the brightest hot Jupiters known, it has no high quality 
phase curve measurements published, and the lack of a wavelength dependent radiative transfer 
routine in our numerical code precludes any detailed comparison with available multi-wavelength 
secondary ec l ipse m easurements. If taken at face value, the observations of HD 209458b by 
Cowan et al.1 (|2007l ) would indicate an infrared flux ratio (F min / F m3iX ) of less than 0.15 — 1.5%; 
however, this is calculated from several discrete obs ervations, a much less reliable way to measure a 
phase curve than cont inuous observation (compare lHarrington et al.1 120061 : ICrossfield et al.l |201Cj) . 
Crossfield et al.l (|2012l ) reported that an existing continuous phase curve observation of HD 209458b 
at 24/im is too plagued by instrumental effects to provide a reliable measurement. 

Fortu nately, there ar e publ ished phase curve results for HD 189733b at a fairly high level of 



Knutson et al.l (120121 ) report multi-wavelength conti nuous phase curve observa tions of 



precision. 

HD 189733b, from 3.6 to 24/j,m (including previous results from iKnutson et~all bool l2009h . since 
hot Jupiters are not blackbodies, but in fact have strong molecular bands in the infrared, there are 
significant differences between the phase curves at each wavelength. The ratio between minimum 
and maximum flux ranges from F m i n /F ma _ x = 15 — 74% and the range for the offset of the hot spot 
from the substellar point is 9p max = 20 - 37°. The work bv lAeol et al.l (|2010h finds values that fall 
within this range and also constrains the temporal variation in 8/im day side emission to be less 
than 2.7% (at 68% confidence). 

Although our model results can only roughly agree with measurements, due our lack of multi- 
wavelength radiative transfer, it is helpful for us to compare the precision of these measurements to 
a mount by which we pre dict that magnetic effects should be able to influence observable properties. 
In IKnutson et al.l (|2012l ) they are able to achieve precisions as good as ±5% for F m \ n / F max and 
±3° for #F max , for observations of HD 189733b. We do not expect magnetic effects to strongly 
influence the circulation on this planet. However, the signature of magnetic effects in the atmosphere 
of HD 209458b is at level comparable to, or even greater than, the precision of the measurements 
for HD 189733b. In other words, if that same precision could be achieved for HD 209458b, the 
magnetic effects would be observable, to the degree that they would need to be included in any 
models used to interpret the observations. 

In addition to the information available from an orbital phase curve, the method of eclipse 
mapping can be used to reconstruct two-dimensional maps of the emission from the day side of 
a planet, by carefully measuring the shape of the light curve as the planet passes into and out of 
secondary eclipse. Recently two gr oups have used this method to create m aps of the 8/xm emission 



from the day side of HD 189733b (|de Wit et al 



2012 



Majeau et al.l 12012 ) . for the first time pro- 



viding information about the latitudinal profile of emission. Both groups confirmed a longitudinal 
offset of the hot spot that is consistent with the phase curve measurement s, but their s e parate 



analysis approaches lead to a difference in their latitudinal profile results, de Wit et al 



found that the hot spot was shifted north of the equator by 17 ± 10°, while iMajeau et al 



(2012) 



20121 ) 



combined the eclipse mapping and phase curve information to arrive at a value consistent with 
no latitudinal shift away from the equator (3.1 ± 9.4° N). A northward shift of the hottest region 
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Table 2. Observable properties of each model 



Model -Fmin/^max #Fmax £ 



HD 209458b, with: 



B 


= G 


17% 


12° 





D 


= 3 G, = 100 K 


13% 


0° 


0.7% 


D 


= 3 G 


13% 


3° 


0.7% 


B 


= 10 G 


12% 


2° 


4% 


B 


= 30 G a 


n/a 


n/a 


10% 



HD 189733b, with: 



B 


= G 


21% 


17° 





B 


= 3 G 


21% 


17° 


0.03% 


B 


= 10 G 


21% 


17° 


0.3% 


B 


= 30 G 


21% 


16° 


2% 


B 


= 30, metallicity x 3 


20% 


16° 


3% 



a This is from the last snapshot before the model crashed 
at 271 P ro t and so is not a robust result, but we include the 
ohmic heating for comparison with other models. 

Note. ■ - The observable properties of each model are: 
the flux ratio between minimum and maximum emission 
(-Pmin/-fmax) 5 the longitude of maximum emission (#Fmax)> 
and the total ohmic heating as a fraction of the stellar input 
(e). In the calculation of e, the total luminosity of incident 
stellar irradiation is 3.3 x 10 22 W for HD 209458b and 1.5 x 
10 22 W for HD 189733b. All quantities are calculated from 
a snapshot of the atmosphere at 2000 P ro t • 
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of the atmosphere would be surprising, given that most m odels of hot Jupiter circulation predict 
strong hemispheric symmetry for the temperature s tructure ( Showman et alj|2009 ; Heng et al. 



Dobbs-Dixon et al.ll2012l ; iRauscher Menoull2012i ). Models that have broken north-south symme- 
try are also those that exhibi t temporal variability (|Cho et al.l 120031 . 120081 ; lLangton &: Laughlin 



most variable rtRauscher et al 



20071 ; iDobbs-Dixon et al.ll201C ) and those able to produce the largest latitu dinal shifts are al so the 



20071 ) and would exceed the limit placed by lAgol et al.1 ()2010l ) . 



As we have already discussed, magnetic effects are unlikely to have any significant influence on 
the atmosphere of HD 189733b, but it is worth considering whether they could in general produce 
observable latitudinal shifts of the hot spot. Upon inspection of several snapshots from each of our 
models, we find strong hemispheric symmetry for all models of HD 189733b, as expected. However, 
the B = 3 and 10 G models of HD 209458b do show variation in the latitude of the hot spot, with 
offsets of up to 10 — 20° N (or S). However, those same models lack any significant offset of the 
hottest region from the substellar longitude. We find that magnetic effects cannot simultaneously 
account for significant offsets away from the substellar point in both latitude and longitude. 



3.6. Global ohmic heating rates 



For many years now it has been recognized that the radii of hot Jupiters should exceed that of 
Jupiter because the intense stellar irradiation will slow their evolutionary cooling and contraction. 
However, many transiting planets still have radii larger than expected, implying an additional 
source of heating that keeps so me planets bloated. S everal mechanisms have been proposed to 
explain these inflated radii (see iFortney et al.l l2010al . for a short review) , but so far none are 
able to explain the full set of observed planets. Ohmic heating is one of the candidates most 
recently proposed to provide the extra source of heating, by tapping into the kinetic energy of the 
atmospheric winds (generated by stellar heating) a nd transporting it deeper into the planet throug h 
the penetration and dissipation of induced currents (jBatygin fc Stevensonhoirt l iPerna et al I hoiobh . 
The success of ohmic heating in producing inflated radii depends on how much power it can generate 
and at what depth, with deeper heating able to have a stronger effect on the planet's evolution 
dGuillot & Showman|[20o3 l. 



One of the difficulties in modeling the evolutionary impact of ohmic heating is that it couples 
regions of the atmosphere with short and long timescales (high in the atmosphere and the deep 
interior). Another is that the one-dimensional models necessary for evolutionary calculations must 
use si mple analytic fo r ms to i nclude the inher e ntly t hree- dimensional nature o f atmo spheric circu- 
lation. iBatvgin et al.l (|201ll ). IWu Lithwick! (|2012l ). and iHuang Cummina (|2012l ) all calculate 
evolutionary models for hot Jupiters that include ohmic heating. They each employ different as- 
sumptions, forms for the wind and temperature profiles, and prescriptions for the effects of feedback 
(such as magnetic drag slowing the winds, and heating of the upper atmosphere), and arrive at 
different conclusions regarding the effectiveness of ohmic heating as a mechanism for explaining the 
observed distribution of hot Jupiter radii. 
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3.6.1. Our calculated heating rates 

Using our models we can calculate the ohmic heating rates directly from the wind and tem- 
perature structure of the atmosphere, consistently including the influences of magnetic drag and 
heating on that atmospheric structure. By continuously updating the local resistivities throughout 
the atmosphere, our model is able to predict the ohmic heating without any use of a prescribed 
form for the drag or wind profile. In Table [2] we report the total ohmic heating in each model, 
summed over the entire domain of the simulations, from 1 mbar to 100 bar. In Figure [TTl we plot 
global profiles of temperature, wind speed, and ohmic heating as a function of pressure, for several 
of our models. These profiles both demonstrate the interplay between various effects and can be 
used as a guide or boundary condition for one-dimensional evolutionary models. 

As we saw in Section [3.31 the cooler temperatures throughout the atmosphere of HD 189733b 
result in weaker magnetic effects, and there is little difference between the non-magnetic and mag- 
netic models of this planetQ Although the HD 189733b B = 3 G model has faster winds than 
the HD 209458b B = 3G model, its lower atmospheric temperatures result in ohmic heating rates 
about two orders of magnitude less than those for HD 209458b. However, when we frame the total 
ohmic heating rates as efficiencies (relative to the stellar input, e, reported in Table [2]), we find that 
some of the models of HD 189733b have heating efficiencies that are comparable to those of the 
HD 209458b models. In particular, the models of HD 189733b with B = 30 G have efficiencies sim- 
ilar to that of the B = 10 G model of HD 209458b. Nevertheless, evolutionary models that include 
ohmic heating find that, for the same heating efficiency, planets with lower effective temperatures 
have smaller ra dii and Jupiter-mas s planets with T e g < 1400 K experience no significant radius 



inflation at all (jBatygin et al.ll201ll ). This is consistent with the observed radius of HD 189733b, 



which is not larger than predicted by evolutionary models with no extra source of heating. 

The global profiles for our HD 209458b models show that most of the planet's ohmic heating 
occurs high in the atmosphere, where the winds are fast. From previous sections we also know that 
most of this heating occurs on the hot day side (e.g., Figure [2]), where the magnetic timescales are 
the shortest. This is an important point because it means that the globally averaged temperature 
profile is not representative of the temperatures at which the heating occurs. If one-dimensional 
models were to use a globally averaged temperature profile to calculate ohmic heating, they would 
overestimate the atmospheric resistivities and underestimate the heating rates. 

The B = 3 and 10 G models have very similar temperature and wind profiles; the factor of ~6 
increase in ohmic heating for the B = 10 G model (integrated over the entire atmosphere) is mainly 
due to the increased magnetic field strength. Although magnetic effects do not significantly change 
the global temperature profiles, compared to the B = G model, they do strongly suppress wind 



7 We only show the B — and 3 G models in Figure 1111 but the temperature profiles for models up to B — 30 G 
are nearly identical, and the wind speeds in the B = 30 G model are decreased by only a few hundred meters per 
second. 
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Fig. 11. — Profiles of temperature, zonal wind speed, and ohmic heating as a function of pressure, 
for several of our models. The temperature profile is the arithmetic mean at each pressure level 
and the wind profile is the root-mean-square. The ohmic heating profile gives the total heating 
rate at each pressure level (see Appendix A for the integration); by adding the contribution from 
each layer (each point) one can integrate over some region of interest. 

speeds, by ~1 km s _1 . The greatest wind reduction occurs for pressures around 1 bar, although 
the winds are slowed as deep as ~20 bar. The strong decrease in winds leads to a dip in ohmic 
heating at ~1 bar, while the heating below this level depends on the temperature structure of the 
deep atmosphere. 
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The ohmic heating profiles for the T nt = 100 K and T; nt = 500 K versions of our B = 3 G 
HD 209458b model rapidly diverge for pressures greater than 1 bar. The heating rate in the 
T; nt = 100 K model drops by 2 orders of magnitude from 1 to 100 bar, while the heating rate in 
the T; nt = 500 K model increases by an order of magnitude over the same pressure range. The 
winds in the T; n t = 100 K model are slower than in the Tj nt = 500 K model, 100 m s _1 instead 
of 200 m s _1 , but the drastic reduction in ohmic heating is much more strongly dependent on the 
huge temperature difference. In Figure [12] we plot analytic, globally averaged temperature-pressure 
profiles for HD 209458b models with T nt from 100 to 500 K, as well as the corr esponding resistivity 
profiles. The temperature profiles use the analytic formalism of iGuillotl (|2010l ). appropriate for our 
models, and the resistivities are calculated from the temperature profiles using Equation [2j Not 
only is there a difference of three orders of magnitude between the resistivities of the T nt = 100 K 
and 500 K models at 100 bar (the lower boundary of our simulations), but these profiles also 
show opposite behavior: from 10 to 100 bar (and deeper) the resistivities of the T nt = 100 K are 
increasing, while the resistivities of the T n t = 500 K model quickly decrease. This difference leads 
to the disparate behavior of the ohmic heating profiles at depth. The current value of T; n t for 
HD 209458b (or any extrasolar planet) is unknown, meaning that the huge difference between 
these two models — a change in the ohmic heating at depth of at least two orders of magnitude — is 
a quantitative demonstration of our ignorance. 
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Fig. 12. — Left: the analytic globally averaged temperature profiles for models of HD 209458b with 
a range of internal heat fluxes, characterized by values for T nt ranging from 100 to 500 K. Right: 
the corresponding profiles of electrical resistivity {rj) , calculated as per Equation [2j 



The internal heat flux of a planet should decrease as it ages and cools, so in some sense the 
comparison between these two models could be seen as between earlier (T nt = 500 K) and later 
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(Tlnt = 100 K) stages of a planet's evolution, and we would therefore expect to find more ohmic 
heating at depth ear lier in a planet' s lifeti me. Note that our res u lts do not necessarily contradict the 
assumption, used in iBatvgin et all (|201lh and lWu Sz Lithwickl (|2012l ). that the magnetic efficiency 
is constant throughout a planet's evolutionjf] Both studies constrain the total heating throughout 
the atmosphere to be constant, but allow for variation of the depth at which heating occurs. In 
fact, due to the slightly higher temperatures and wind speeds in the upper atmosphere of our 
Tint = 100 K model, it has more ohmic heating at altitude than the Ti nt = 500 K model. This 
compensates for the decrease in heating at deep pressures and results in these two models having 
the same total rate of ohmic heating (see Table [2]). Since we only have this single example of 
identical models with different T{ nt , we cannot comment on whether this is somehow a fundamental 
property, or just coincidence. 



3.6.2. Our findings on ohmic heating as a cause for HD 209458b's inflated radius 



Finally, we consider whether the amount of ohmic heating in any of our HD 209458b m odels 
is sufficient to explain this planet's bloated radius. According to iGuillot &: Showman! (120021 ). the 
inflated radius of HD 209458b requires roughly 10% of the stellar insolation to be deposited as 
heating around the 5 bar level, e = 1% if around 20 bar, or e ~ 0.08% if in the adiabatic interior, 
which begins at ~160 bar in their model. On the other hand, the results of lBatygin Stevenson 
(|2010l ) are able to reproduce the radius of HD 209458b with 4 x 10 18 W of ohmic heating (or 
e ~ 0.01%) in the interior, which in their model begins at ~90 bar, for the case of solar metallicity, 
no core, and Tj so = 1700 K. Note that the amount of heating required depends not just on the 
pressure at which the heating is deposited, but also on the details of the particular model, including 
the metallicity of the atmosphere, the presence or absence of a core, the location of the radiative- 



connective boundary, and how all these vary througho ut the planet's evolution (see iBatygin et al 



2011 



Huang fe Cummingj 120121 : IWu &: Lithwickl l2012l ) . This motivates future work to develop a 



more consistent connection between ohmic heating results from evolutionary and circulation models. 



Our B = 3 G model cannot fulfill the IGuillot fc Showman! (|2002l ) requirement of e = 1% at 
20 bar, nor can the B = 10 G model, whose ohmic heating, integrated from 10 to 100 bar, only 
amounts to ~0.2% of the stellar insolation. Although the total ohmic heating in this model is 4% 
of the stellar heating, most of it is deposited at altitude and cannot effect the planet's evolution. 
Due to the use of constant absorption coefficients in our radiative transfer scheme, the temperature 
profiles in our model will never become adiabatic (see Appendix B of RM12) and so we have no 
radiative-convective boundary. If we consider the ohmic power in the deepest level of our models, 
which spans ~ 70 — 100 bar, then both our B = 10 G and B = 3 (T| n t = 500 K) models easily fulfill 



8 While lHuang fc Cumminel |2012h do not assume a constant efficiency, their models have ohmic heating at depth 
that is not monotonic with the inner entropy and it is difficult to make a direct comparison between their models 
and ours. 
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the iBatvgin fc Stevenson! (|2010l ) requirement of 4 x 10 18 W at ~90 bar. However, the B = 3 G 
model with Tj nt = 100 K is orders of magnitude away from reaching that heating threshold, again 
demonstrating the importance of the deep thermal structure on the amount of ohmic heating. 



We can also employ the scalings from IWu &: Lithwickl (120121) t o estimate the ohmic power 
dissipated below our models' bottom boundary. IWu fe Lithwickl (120121 ) argue that the conservation 
of currents implies simple scalings between the radial and meridional components of atmospheric 
currents and those in the interior. They make the case, based on conservation of total current, 
that at an order of magnitude level the total heating in the interior must be ~ z w i n d/R p smaller 
than in the "weather layer" (which has a vertical extent z w i n d)- The average thickness of our 
modeled atmosphere, from 100 bar to 1 mbar, is ~ 2 x 10 6 m for all of our HD 209458b models. 
From z w [ n d/R p = 0.02 and the total ohmic heating rates in each of our models we find that the 
heating below our bottom boundary should be ~ 6 x 10 18 W for B = 3 G and ~ 3 x 10 19 W 
for = 10 G. Both estimates are above the 4 x 10 18 W threshold from IBatvgin &: Stevenson 



2010) and the B = 10 G model matches the interior heating constraint from lGuillot Showman 



20021 ). (I t is worth noting that our he ating rates are not too dissimilar from the solar composition 
models of IBatvgin Stevenson! (|2010l ) with Tj so = 1700 K, with a few percent for the total ohmic 
heating effi ciency but only about a few times 10 19 W dissipated deeper than 10 bar.) From the 
scalings of IWu Lithwickl (|2012l ). we would calculate identical interior dissipation rates for the 
Tint = 100 K and T- m i = 500 K models with B = 3G, although we have shown that ohmic heating 
at depth depends very strongly on the deep atmosphere profile. This is clearly a limitation of their 
idealized model, which assumes a constant conductivity from the upper to deep atmosphere, with 
the conductivity only changing at the transition to the interior. Again, this is an issue that would 
benefit from more complete models that couple the atmosphere and interior. 



The estimate from the scalings in lWu Lithwickl (|2012l ) for our interior heating rates rests on 
radial currents alone, while the formalism we use for the drag only describes meridional currents 
(which are dominant in the thin induction region). In this way the theories are complementary, 
since meridional and radial currents are additive in terms of ohmic dissipation. Therefore, our 
estimates from this scaling argument, together with the rising ohmic power at depth in our models 
(based on meridional currents), both suggest that interior ohmic dissipation can inflate HD 209458b 
for planetary magnetic field strengths of B > 3 — 10 G. 



4. Summary 

We have updated our three-dimensional model of atmospheric circulation to include geomet- 
rically and energetically consistent magnetic drag and ohmic heating. This is the first model to 
directly couple the magnetic effects to the full atmospheric structure. We calculate the electrical 
resistivities directly from local conditions and update these values at each timestep in the simula- 
tion. We apply magnetic drag only to the zonal (east-west) component of the flow and include the 
latitudinal dependence in the strength of the drag. All kinetic energy lost through magnetic drag 
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is consistently returned to the atmosphere as localized ohmic heating. In these ways the magnetic 
effects are strongly coupled with the atmospheric circulation. Here we present results from this 
code for two well known hot Jupiters, HD 189733b and HD 209458b. We test models of these 
planets with planetary magnetic field strengths ranging from to 30 G. 

Due to the ~300 K difference in equilibrium temperature between HD 189733b and HD 290458b, 
magnetic drag and ohmic heating only influence the circulation of the hotter planet, HD 209458b. 
The B = 3G model of HD 209458b is obviously different from the non-magnetic version, but even 
the B = 30 G model of HD 189733b with enhanced metallicity does not show significant differences 
from the B = G model, having only slightly slower wind speeds. The models of HD 189733b 
also do not show any signature of magnetic effects in their observable properties, nor do they pro- 
duce enough ohmic heating to effect the planet's evolution, which is consistent with the observed 
(non-inflated) radius of this planet. 

We find that magnetic effects are able to strongly influence the circulation of HD 209458b; in 
several aspects the models with B > differ from the B = version. In particular, the magnetic 
models: 1) have slower wind speeds, by ~1 km s^ 1 , 2) do not have an equatorial eastward jet 
that circles the globe, 3) have departures from hemispheric symmetry in the temperature and flow 
patterns, and 4) maintain more of a hot-day/cold- night temperature structure, over a deeper range 
of pressures, than the non-magnetic model. These trends gradually become stronger at higher 
magnetic field strengths. We also find that it is the magnetic drag that has a stronger influence on 
the circulation than the (coupled) ohmic heating, both of which act primarily on the day side of the 
planet. Throughout almost all of the atmosphere the local ohmic heating is a very small fraction 
of the radiative heating. However, our B = 30 G model of HD 209458b has localized features 
where the ohmic heating reaches as much as 10% of the local radiative heating; this model became 
numerically unstable and crashed before completion. 

Magnetic drag and heating have a strong enough effect on the circulation of HD 209458b that 
they alter its observable properties. The flux contrast between the day and night sides of the 
planet is greater in the magnetic models; during the planet's orbit the minimum flux emitted is 
only 12 — 13% of the maximum flux, compared to a ratio of 17% for the non-magnetic model. 
We also find that the brightest region of the atmosphere remains well aligned with the substellar 
longitude, compared to a 12° eastward shift in the B = G model. Although there is no shift in 
longitude in the magnetic models, we do find that the latitude of the brightest region varies in time 
and can shift away from the substellar point by as much as 10 — 20°. The differences between the 
magnetic and non-magnetic models of HD 209458b are at a level such that they could be measured, 
meaning that models used to interpret observations of this planet should include magnetic effects. 

Finally, we compare the ohmic heating profiles from our models of HD 209458b to predictions 
from evolutionary models to determine whether we find sufficient heating at depth to explain the 
inflated radius of this planet. Most of the ohmic heating in our models occurs high in the atmosphere 
and cannot prevent the planet's standard cooling and contraction. However, the heating in the 



-32- 



deepest layers of our B = 3 and 10 G models could fulfill the requirement for inflation set by 



Batvgin fc Stevenson! (|2010i ). although only if the internal heat flux is high enough for a hot deep 
atmosphere; we found heating rates at depth to be two orders of magnitude higher in the B = 3G 
model with T[ n t = 500 K than the one with T[ Qt = 100 K. The models with T; n t = 500 K also 
have heating profiles that are increasing with pressure from 10 to 100 bar, even though the winds 
are becoming s l ower, with speeds on the order of 100 m s _1 . We use scaling arguments from 
Wu &: Lithwickl ([2012J) to estimate the ohmic heating in the planet's interio r, below our models' 
bottom bo undary. These estima t es m eet or exceed the requirements from iGuillot Showman 
(|2002l ) and iBatvgin fc Stevenson! (|2010l ) for ohmic heating in the adiabatic interior. Both these 
estimates, and the rising ohmic power with pressure that we find in our models, suggest that interior 
ohmic dissipation can inflate the radius of HD 209458b for planetary magnetic field strengths at or 
greater than B = 3 — 10 G. 

In order to model the complex interaction between magnetic effects and the atmospheric cir- 
culation, we made several simplifying assumptions. We expect these choices to be appropriate to 
first order, but as always, we would benefit from a fuller theory, in which these assumptions could 
be relaxed. In particular, our calculation of the magnetic drag and ohmic heating is b ased on a 



form alism that assumes axisymmetry in the flow structure and atmospheric resistivities ( Liu et al 



20081 ). neither of which is realized in hot Jupiter atmospheres, especially at low pressures. We 
have also assumed the simplest magnetic geometry, with an aligned dipole field, while in reality the 
planet's magnetic axis could be misaligned from its axis of rotation, or the field could be multipole 
or uneven in other ways. Varying the magnetic geometry could have strange impacts on the circu- 
lation. In addition, our magnetic drag is only applied to the zonal flow, although the meridional 
flow may also begin to experience drag as it nears the poles and the field becomes more radial, 
but we lack the formalism to include this effect. However, even with these simplifications, we have 
established that the inclusion of magnetic effects results in a greater richness of possibilities for 
atmospheric dynamics on the hottest exoplanets. 



This work was performed in part under contract with the California Institute of Technology 
(Caltech) funded by NASA through the Sagan Fellowship Program. ER thanks the Physics Depart- 
ment at Virginia Tech for kindly hosting her during the writing of this paper. KM was supported 
by the NASA grant PATM NNX11AD65G. 



A. An explanation of the different forms of ohmic heating shown in plots 

In Figures [21 HI and [6] we plot the specific ohmic heating rate in units of W kg . These rates 
are calculated from the local atmospheric structure, following Equations [H through [6l 

In Figure [7] we plot the ohmic heating rate as a function of latitude for various models. Here 
we have calculated the heating rates everywhere in the model (as per above) and then summed 
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over longitude and pressure at each latitude: 

Q(<£)/(#) = ! {u 2 /T mag ) dm/dcj) = ( (n 2 /r mag ) (l/g)dP fljcos(0) d(9 (Al) 
J Jp,e 

(where we have taken advantage of hydrostatic balance to convert dm = pdzdA to pressure coor- 
dinates). The resulting heating rates are in units of watts per radian. 

In Figure HU we plot the total ohmic heating at each pressure level, in units of watts: 
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where AP is the thickness of the level in units of pressure. 
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